Advanced Numerical Methods II

Vadim Sokolov Northern Illinois University sokolov@math.niu.edu

Wolfram Research  Mathematica Technology Conference, Champaign, IL- 10/2005  Introduction

The package works with a finite-dimensional time-invariant linear continous-time dynamical system with the following representation:

\begin{matrix}ẋ(t) = Ax(t) + Bu(t), t ≥ 0 \\ y(t) = Cx(t) + Du(t) \\ x(0) = x_{0}, \end{matrix}  (1)

where A ∈ ℝ^{n×n}, B ∈ ℝ^{n×m}, C ∈ ℝ^{r×n}, D ∈ ℝ^{r×m}, n is the order of (1). Assumed, that m << n

Wolfram Research  Mathematica Technology Conference, Champaign, IL- 10/2005  Capabilities

• Solution of continuous-time algebraic Lyapunov equations, i.e equation of the type

AX + X^{T} A + BB^{T} = 0, A ∈ ℝ^{n×n}, B ∈ ℝ^{n×m}  (2)

• Model Reduction The idea of model reduction is to construct a reduced-order model from the original full-order model

\begin{matrix}ẋ(t) = Ãx(t) + B̃x(t), t ≥ 0 \\ y(t) = C̃x(t) + D̃u(t) \\ x(0) = x_{0}, \end{matrix}  (3)

such that the reduced-order model is close, in some sense, to the original model

• Solution of Algebraic Riccati Equation:

XA + A^{T} X − XBR^{−1}B^{T} X + C^{T} QC = 0  (4)

Wolfram Research  Mathematica Technology Conference, Champaign, IL- 10/2005  Motivation

With availability of LAPACK it is possible to handle control problems of order 1000, using classical methods. Classical methods implemented in CSP-ANM destroy the sparsity, inherited in most large practical problems Sparsity is lost by the use of canonical forms such as, triangular, Hessenberg and real-Schur, which are obtained by using Gaussian eliminations, Householder and Givens transformations.

Wolfram Research Mathematica Technology Conference, Champaign, IL- 10/2005 Solution of Lyapunov Equation

• Krylov subspace method INPUT: A, B, m OUTPUT: X̃ - approximate solution to Lyapunov Equation

1: {Q,R} = QRDecomposition[B] 2: run block – arnoldi algorithm for matrix A, with starting block – vector Q, obtain

H_{m}, V_{m}

3: solve projected Lyapunov equation

H_{m}G_{m} + G_{m}H_{m}^{T} + ( \begin{matrix}R \\ 0 \end{matrix} ) ( R^{T} 0 ) = 0

4: Return X̃ = V_{m}G_{m}V_{m}^{T} • Low rank Cholesky factor ADI iteration for the Lyapunov equation INPUT: A, B OUTPUT: X̃ - approximate solution to Lyapunov Equation 1: compute iteration parameters, (p_{1}, …, p_{J})

2: Define P_{i} = (^{\begin{matrix}√ & \overline{−2p_{i+1}} \\ √ & \overline{−2p_{i}} \end{matrix}})[I − (p_{i+1} + p_{i})(A + p_{i+1}I)^{−1}]

3: z_{1} = (\sqrt{−2p_{1}})(A + p_{1}I)^{−1}B 4: Z_{1} = [z_{1}] 5: for j = 1, …, J do 6: z_{j} = P_{j−1}z_{j−1} 7: Z_{j} = [Z_{j−1} z_{j}]

8: end for

If we take X_{0} to be zero matrix then X_{j} is symmetric.

Wolfram Research  Mathematica Technology Conference, Champaign, IL- 10/2005  Matrix Operations

• Solution of a linear system with with shifted matrix A

• Matrix - vector multiplication

Wolfram Research  Mathematica Technology Conference, Champaign, IL- 10/2005  When to use?

• Size of the problem is > 500

• Matrix A has to be stable, i.e all eigenvalues lie in the left half-plane

• Number of inputs and outputs is considerably less then the number of states, i.e

m << n

• Matrix A has a structure, i.e computation of matrix-vector product and solution of linear system with matrix A are relatively cheap operations

Wolfram Research  Mathematica Technology Conference, Champaign, IL- 10/2005  Test Example

Consider the following PDE:

ẋ = \frac{∂^{2}x}{∂z^{2}} + \frac{∂^{2}x}{∂v^{2}} + 20\frac{∂x}{∂z} − 180x + F (z, v, t)  (5)

x = x(z, v, t), assume that F (z, v, t) has the form

F (z, v, t) = f(z, v)g(t)

If we discretize the region with n interior points in the y direction and m interior points in the v direction, then the resulting matrix problem is of the form

ẋ = Ax + Bu,

Where A is an nm by nm matrix. We put output matrix C to be equal to B, for the reason of simplicity.

Wolfram Research  Mathematica Technology Conference, Champaign, IL- 10/2005  Future Plans

• Parameter - free Krylov projection method for Lyapunov Equation

• Krylov Projection Methods for model reduction

• Krylov Projection Methods for Riccati Equation

• Partial Pole assignment

• Documentation